#DESCRIPTION
#Defining slice parameters for a given slope material and geometry, and slip sirface. For use in limit equilibrium analysis

#1. libraries and modules
import numpy as np
import pandas as pd
import math

#2. math/trigo functions
tan=math.tan
atan=math.atan
rad=math.radians

#3. determining coordinates of top and bottom edges of slices
def slice_pts(H,B,h,k,R,x1,x2,num_slices):
    #INPUTS
    #H,B:           floats, slope height (H) and gradient (B)
    #num_slices:    integer; number of slices to be used in the limit equilibrium analysis
    #h,k,R:         floats; circular slip surface center (h,k) and radius (R)
    #x1, x2:        floats; x-coordiante of slip surface endpoints
    
    x=np.linspace(x1,x2,num_slices+1)
    yt=np.where(x>0,
                np.where(x>H/tan(rad(B)),
                         H*np.ones(len(x)),
                         x*tan(rad(B))),
                np.zeros(len(x)))

    yb=np.where((h-x)**2==R**2,
                np.ones(len(x))*k,
                (k-np.sqrt(R**2-(h-x)**2)))

    return x,yt,yb

#3. defining slice parameters independent of LE analysis
def Slice_params(H,B,num_slices,h,k,R,x1,x2,uw,ru,kh, full_output_Slices):
    #INPUTS
    #H,B:           floats, slope height (H) and gradient (B)
    #num_slices:    integer; number of slices to be used in the limit equilibrium analysis
    #h,k,R:         floats; circular slip surface center (h,k) and radius (R)
    #x1, x2:        floats; x-coordiante of slip surface endpoints
    #uw:            float; unit weight of slope material
    #ru:            float; pore water pressure coefficient
    #kh:            float; horizontal seismic coefficient

    

    #1. determining coordinates of top and bottom edges of slices
    x,yt,yb=slice_pts(H,B,h,k,R,x1,x2,num_slices)
    
    #2. computing width of each slice
    width=x[1]-x[0]

    #3. other parameters
    height=[]               #height of slice, m
    alpha=[]                #inclination of slice base, radians
    v_seis=[]               #vertical distance of horizontal seismic force from circle center, m
    pwp=[]                  #pore water pressure at slice base, MPa
    z=[]                    #height of water table above slice base
    v_tcF=[]                #vertical distance of tension crack force from circle center, m
    uww=9.81*10**-3         #unit weight of water, MN/m3
    for i in range(len(x)-1):
        height.append(0.5*((yt[i]+yt[i+1]) - (yb[i]+yb[i+1])))
        alpha.append((atan((yb[i+1]-yb[i])/width)))
        v_seis.append(k-0.5*(yb[i+1]+yb[i]+height[i]))
        pwp.append(height[i]*uw*ru)
        z.append(pwp[i]/uww)
        v_tcF.append(k-(yb[i]+(z[i]/3.)))
    height=np.asarray(height)
    alpha=np.asarray(alpha)
    v_seis=np.asarray(v_seis)
    pwp=np.asarray(pwp)
    z=np.asarray(z)
    v_tcF=np.asarray(v_tcF)

    bl=width/np.cos(alpha)          #length of slice base, m
    weight=height*width*uw          #weight of slice, MN
    SF=weight*np.sin(alpha)         #shearing force due to slice weight, MN
    seisF=kh*weight*v_seis/R        #seismic force acting on slice, MN
    tcF=(0.5*uww*z**2)*v_tcF/R      #tension crack force acting on slice, MN

    if full_output_Slices==1:
        Slices=pd.DataFrame(data=None,index=range(1,num_slices+1))
        Slices['height']=height
        Slices['width']=width*np.ones(len(height))
        Slices['alpha']=np.degrees(alpha)
        Slices['baseL']=bl
        Slices['weight']=weight

        Slices['SF']=SF
        Slices['seisF']=seisF
        Slices['pwp']=pwp
        Slices['z']=z
        Slices['TcF']=tcF
    else:
        Slices=[]

    return weight,alpha,bl,pwp,SF,seisF,tcF, Slices
